# Nick Dietrich
# Explaining Support for International Action Against Human Rights Abusers
# Replication Materials
# dietrich.nicholas@gmail.com

# clear memory
rm(list=ls())

# set working directory, log results, and set seed
setwd("~/HR Support Replication") # path to folder containing dietrich_replication.csv
#sink("dietrich_replication_log.txt")
set.seed(666)

# read in data
dat <- read.csv("dietrich_replication.csv")

### treatment variables, coded from least precise to most precise
# num: precision in number of victims
# loc: precision in location
# time: promptness of reporting (year, month, day)
### outcome variables, coded from least to most support
# denounce: support for denouncing perpetrator
# sanction: support for economic sanctions
# interven: support for military intervention
### control variables
# country: country of residence of respondent
# attn: respondent passed attention check (binary)
# liberal: respondent's self-reported ideology, from most conservative (-3) to most liberal (3)

###
### Table 2: Linear regression models, with and without controls, for each outcome
###
ols.denounce1 <- lm(denounce~as.factor(num)+as.factor(time)+as.factor(loc),data=dat)
summary(ols.denounce1)
ols.denounce2 <- lm(denounce~as.factor(num)+as.factor(time)+as.factor(loc)+
                        as.factor(country=="US")+liberal,data=dat)
summary(ols.denounce2)
ols.sanction1 <- lm(sanction~as.factor(num)+as.factor(time)+as.factor(loc),data=dat)
summary(ols.sanction1)
ols.sanction2 <- lm(sanction~as.factor(num)+as.factor(time)+as.factor(loc)+
                        as.factor(country=="US")+liberal,data=dat)
summary(ols.sanction2)
ols.interven1 <- lm(interven~as.factor(num)+as.factor(time)+as.factor(loc),data=dat)
summary(ols.interven1)
ols.interven2 <- lm(interven~as.factor(num)+as.factor(time)+as.factor(loc)+
                        as.factor(country=="US")+liberal,data=dat)
summary(ols.interven2)

###
### A brief sensitivity analysis: how large would our sample need to be to find significant effects?
###
denounce.robust1 <- lm(denounce~as.factor(num)+as.factor(time)+as.factor(loc),data=rbind(dat,dat,dat,dat))
summary(denounce.robust1)
denounce.robust2 <- lm(denounce~as.factor(num)+as.factor(time)+as.factor(loc)+
                      as.factor(country=="US")+liberal,data=rbind(dat,dat,dat,dat))
summary(denounce.robust2)
sanction.robust1 <- lm(sanction~as.factor(num)+as.factor(time)+as.factor(loc),data=rbind(dat,dat,dat,dat))
summary(sanction.robust1)
sanction.robust2 <- lm(sanction~as.factor(num)+as.factor(time)+as.factor(loc)+
                      as.factor(country=="US")+liberal,data=rbind(dat,dat,dat,dat))
summary(sanction.robust2)
interven.robust1 <- lm(interven~as.factor(num)+as.factor(time)+as.factor(loc),data=rbind(dat,dat,dat,dat))
summary(interven.robust1)
interven.robust2 <- lm(interven~as.factor(num)+as.factor(time)+as.factor(loc)+
                      as.factor(country=="US")+liberal,data=rbind(dat,dat,dat,dat))
summary(interven.robust2)

# answer: more than four times as large, except for the wrong-direction effect of the number treatment.

###
### Table 2: OLS with confidence intervals
###
library(stargazer)
out <- "tab_ci.tex" # ordinal
star <- stargazer(ols.denounce1,ols.denounce2,ols.sanction1,ols.sanction2,ols.interven1,ols.interven2,
                  align=T,digits=2,star.cutoffs=.05,keep.stat=c("n","AIC"),ci=T,
                  column.labels=c("Denouncing at U.N.","Imposing Sanctions","Sending Troops"),
                  column.separate=c(2,2,2),model.numbers=F,
                  notes="Support for each outcome ranges from zero, indicating strong opposition, to three, indicating strong support. 95% confidence intervals are in parentheses. *p<0.05",
                  notes.append=F,notes.align="l",
                  font.size = "small",column.sep.width = "0pt",
                  title="Support for Action Against Perpetrator, Linear Regression",
                  label="tab-ci",no.space=T,ord.intercepts=F,
                  covariate.labels=c("Exact Number","One Month Ago","Yesterday",
                                     "Exact Location","U.S. Resident","Liberal Ideology"),
                  dep.var.labels=c("Support for","Support for","Support for"))
note.latex <- "\\multicolumn{7}{l} {\\parbox[t]{\\textwidth}{ \\textit{Notes:} Support for each outcome ranges from zero, indicating strong opposition, to three, indicating strong support. 95\\% confidence intervals are in parentheses. *p<0.05.}} \\\\"
star[grepl("Note",star)] <- note.latex
cat(star, sep = '\n', file = out)

###
### Table C: Sensitivity analysis (appendix)
###

### write output as table
library(stargazer)
out <- "tab_sensitivity.tex" # ordinal
star <- stargazer(denounce.robust1,denounce.robust2,sanction.robust1,sanction.robust2,interven.robust1,interven.robust2,
                  align=T,digits=2,star.cutoffs=.05,keep.stat=c("n","AIC"),
                  column.labels=c("Denouncing at U.N.","Imposing Sanctions","Sending Troops"),
                  column.separate=c(2,2,2),model.numbers=F,
                  notes="Support for each outcome ranges from zero, indicating strong opposition, to three, indicating strong support. Standard errors are in parentheses. *p<0.05",
                  notes.append=F,notes.align="l",
                  title="Linear Regression Sensitivity Analysis with Simulated Large Sample",
                  label="tab-sensitivity",no.space=T,ord.intercepts=F,
                  covariate.labels=c("Exact Number","One Month Ago","Yesterday",
                                     "Exact Location","U.S. Resident","Liberal Ideology"),
                  dep.var.labels=c("Support for","Support for","Support for"))
note.latex <- "\\multicolumn{7}{l} {\\parbox[t]{\\textwidth}{ \\textit{Notes:} Samples are simulated by multiplying respondents from the original sample. Results demonstrate the effect of increasing N on the statistical significance of coefficients, assuming that the original point estiamtes remain unchanged. Support for each outcome ranges from zero, indicating strong opposition, to three, indicating strong support. Standard errors are in parentheses. *p<0.05.}} \\\\"
star[grepl("Note",star)] <- note.latex
cat(star, sep = '\n', file = out)


###
### Table A: ordinal logistic regression, with and without controls, for each outcome
###
library(MASS)
ord.denounce1 <- polr(as.factor(denounce)~as.factor(num)+as.factor(time)+as.factor(loc),data=dat,method="logistic")
summary(ord.denounce1)
ord.denounce2 <- polr(as.factor(denounce)~as.factor(num)+as.factor(time)+as.factor(loc)+
                        as.factor(country=="US")+liberal,data=dat,method="logistic")
summary(ord.denounce2)
ord.sanction1 <- polr(as.factor(sanction)~as.factor(num)+as.factor(time)+as.factor(loc),data=dat,method="logistic")
summary(ord.sanction1)
ord.sanction2 <- polr(as.factor(sanction)~as.factor(num)+as.factor(time)+as.factor(loc)+
                        as.factor(country=="US")+liberal,data=dat,method="logistic")
summary(ord.sanction2)
ord.interven1 <- polr(as.factor(interven)~as.factor(num)+as.factor(time)+as.factor(loc),data=dat,method="logistic")
summary(ord.interven1)
ord.interven2 <- polr(as.factor(interven)~as.factor(num)+as.factor(time)+as.factor(loc)+
                        as.factor(country=="US")+liberal,data=dat,method="logistic")
summary(ord.interven2)

### add aic to ordinal models
ord.denounce1$AIC <- AIC(ord.denounce1)
ord.denounce2$AIC <- AIC(ord.denounce2)
ord.sanction1$AIC <- AIC(ord.sanction1)
ord.sanction2$AIC <- AIC(ord.sanction2)
ord.interven1$AIC <- AIC(ord.interven1)
ord.interven2$AIC <- AIC(ord.interven2)

### write output as table
library(stargazer)
out <- "tab_ord.tex" # ordinal
star <- stargazer(ord.denounce1,ord.denounce2,ord.sanction1,ord.sanction2,ord.interven1,ord.interven2,
                  align=T,digits=2,star.cutoffs=.05,keep.stat=c("n","AIC"),
                  column.labels=c("Denouncing at U.N.","Imposing Sanctions","Sending Troops"),
                  column.separate=c(2,2,2),model.numbers=F,
                  notes="Support for each outcome is coded as an ordinal variable that ranges from zero, indicating strong opposition, to three, indicating strong support. Standard errors are in parentheses. *p<0.05",
                  notes.append=F,notes.align="l",
                  title="Support for Action Against Perpetrator, Ordinal Logistic Regression",
                  label="tab-ord",no.space=T,ord.intercepts=F,
                  covariate.labels=c("Exact Number","One Month Ago","Yesterday",
                                     "Exact Location","U.S. Resident","Liberal Ideology"),
                  dep.var.labels=c("Support for","Support for","Support for"))
note.latex <- "\\multicolumn{7}{l} {\\parbox[t]{\\textwidth}{ \\textit{Notes:} Support for each outcome is coded as an ordinal variable that ranges from zero, indicating strong opposition, to three, indicating strong support. Standard errors are in parentheses. *p<0.05.}} \\\\"
star[grepl("Note",star)] <- note.latex
cat(star, sep = '\n', file = out)

###
### Table B: binary logistic regression, with and without controls, for each outcome
###
logit.denounce1 <- glm(denounce.d ~ as.factor(num) + as.factor(time)
                       + as.factor(loc), data=dat, family="binomial")
summary(logit.denounce1)
logit.denounce2 <- glm(denounce.d ~ as.factor(num) + as.factor(time)
                       + as.factor(loc) + as.factor(country=="US")+liberal, 
                       data=dat, family="binomial")
summary(logit.denounce2)
logit.sanction1 <- glm(sanction.d ~ as.factor(num) + as.factor(time)
                       + as.factor(loc), data=dat, family="binomial")
summary(logit.sanction1)
logit.sanction2 <- glm(sanction.d ~ as.factor(num) + as.factor(time)
                       + as.factor(loc) + as.factor(country=="US")+liberal, 
                       data=dat, family="binomial")
summary(logit.sanction2)
logit.interven1 <- glm(interven.d ~ as.factor(num) + as.factor(time)
                       + as.factor(loc), data=dat, family="binomial")
summary(logit.interven1)
logit.interven2 <- glm(interven.d ~ as.factor(num) + as.factor(time)
                       + as.factor(loc) + as.factor(country=="US")+liberal, 
                       data=dat, family="binomial")
summary(logit.interven2)

### write output as table
out <- "tab_logit.tex"
star <- stargazer(logit.denounce1,logit.denounce2,logit.sanction1,logit.sanction2,logit.interven1,logit.interven2,
                  align=T,digits=2,star.cutoffs=.05,keep.stat=c("n","AIC"),
                  column.labels=c("Denouncing at U.N.","Imposing Sanctions","Sending Troops"),
                  column.separate=c(2,2,2),model.numbers=F,
                  notes="Support for each outcome is coded as a binary variable where zero indicates opposition or strong opposition and one indicates support or strong support. Standard errors are in parentheses. *p<0.05",
                  notes.append=F,
                  title="Support for Action Against Perpetrator, Logistic Regression",
                  label="tab-logit",no.space=T,
                  covariate.labels=c("Exact Number","One Month Ago","Yesterday",
                                     "Exact Location","U.S. Resident","Liberal Ideology"),
                  dep.var.labels=c("Support for","Support for","Support for"))
note.latex <- "\\multicolumn{7}{l} {\\parbox[t]{\\textwidth}{ \\textit{Notes:} Support for each outcome is coded as a binary variable where zero indicates opposition or strong opposition and one indicates support or strong support. Standard errors are in parentheses. *p<0.05.}} \\\\"
star[grepl("Note",star)] <- note.latex
cat(star, sep = '\n', file = out)


###
### Figure 3: Mean level of support by treatment group
###
# Calculate bootstrapped CIs
n <- 1000 # number of bootstrapped samples
boot.denounce <- data.frame(matrix(nrow=n,ncol=7)) # 7 treatment levels with 3 outcomes = 21 parameters of interest
boot.sanction <- data.frame(matrix(nrow=n,ncol=7))
boot.interven <- data.frame(matrix(nrow=n,ncol=7))
colnames(boot.denounce) <- c("num_imprecise","num_precise","time_year","time_month","time_now","loc_imprecise","loc_precise")
colnames(boot.sanction) <- c("num_imprecise","num_precise","time_year","time_month","time_now","loc_imprecise","loc_precise")
colnames(boot.interven) <- c("num_imprecise","num_precise","time_year","time_month","time_now","loc_imprecise","loc_precise")

# Loop through, generate bootstrapped samples, and save the means
for(i in 1:n){
  # resample with replacement:
  tmp <- dat[sample(x=1:nrow(dat),size=nrow(dat),replace=T),]
  # store means conditional on treatment conditions
  # support for denouncing perpetrator:
  boot.denounce$num_imprecise[i] <- mean(tmp$denounce[which(tmp$num==0)])
  boot.denounce$num_precise[i] <- mean(tmp$denounce[which(tmp$num==1)])
  boot.denounce$time_year[i] <- mean(tmp$denounce[which(tmp$time==0)])
  boot.denounce$time_month[i] <- mean(tmp$denounce[which(tmp$time==1)])
  boot.denounce$time_now[i] <- mean(tmp$denounce[which(tmp$time==2)])
  boot.denounce$loc_imprecise[i] <- mean(tmp$denounce[which(tmp$loc==0)])
  boot.denounce$loc_precise[i] <- mean(tmp$denounce[which(tmp$loc==1)])
  # support for economic sanctions:
  boot.sanction$num_imprecise[i] <- mean(tmp$sanction[which(tmp$num==0)])
  boot.sanction$num_precise[i] <- mean(tmp$sanction[which(tmp$num==1)])
  boot.sanction$time_year[i] <- mean(tmp$sanction[which(tmp$time==0)])
  boot.sanction$time_month[i] <- mean(tmp$sanction[which(tmp$time==1)])
  boot.sanction$time_now[i] <- mean(tmp$sanction[which(tmp$time==2)])
  boot.sanction$loc_imprecise[i] <- mean(tmp$sanction[which(tmp$loc==0)])
  boot.sanction$loc_precise[i] <- mean(tmp$sanction[which(tmp$loc==1)])
  # support for military intervention
  boot.interven$num_imprecise[i] <- mean(tmp$interven[which(tmp$num==0)])
  boot.interven$num_precise[i] <- mean(tmp$interven[which(tmp$num==1)])
  boot.interven$time_year[i] <- mean(tmp$interven[which(tmp$time==0)])
  boot.interven$time_month[i] <- mean(tmp$interven[which(tmp$time==1)])
  boot.interven$time_now[i] <- mean(tmp$interven[which(tmp$time==2)])
  boot.interven$loc_imprecise[i] <- mean(tmp$interven[which(tmp$loc==0)])
  boot.interven$loc_precise[i] <- mean(tmp$interven[which(tmp$loc==1)])
  # remove temporary dataset
  rm(tmp)
}

# plot the actual sample means with bootstrapped CIs
library(tikzDevice)
tikz("boot_means.tex",height=5.5,width=6.5)
#pdf("boot_means.pdf",height=5.5,width=6.5)
par(mfrow=c(3,1),mar=c(3,9.2,2,1))
### Denounce
plot(c(NA,NA),xlim=c(-.1,3.1),ylim=c(.9,9.1),main="Denounce at UN",axes=F,xlab="",ylab="")
#boxplot(x=dat$denounce,xlim=c(-.1,3.1),ylim=c(.9,3.1),main="Denounce at UN",axes=F,xlab="",ylab="",horizontal=T)
axis(1,at=c(0,1,2,3),labels=c("Strongly Oppose","Oppose","Support","Strongly Support"))
abline(v=quantile(dat$denounce,.25),lty=2,col="grey80") # quantile, all treatment groups
abline(v=quantile(dat$denounce,.75),lty=2,col="grey80") # quantile, all treatment groups
abline(v=mean(dat$denounce),lty=2) # sample mean, all treatment groups
## NUMBER
# imprecise:
lines(x=c(quantile(boot.denounce$num_imprecise,probs=.025),quantile(boot.denounce$num_imprecise,probs=.975)),y=c(8,8))
points(mean(dat$denounce[which(dat$num==0)]),8,pch=19)
axis(2,at=c(8,9),labels=c("Unknown Number","Exact Number"),las=1)
# precise:
lines(x=c(quantile(boot.denounce$num_precise,probs=.025),quantile(boot.denounce$num_precise,probs=.975)),y=c(9,9))
points(mean(dat$denounce[which(dat$num==1)]),9,pch=19)
## TIME
# year
lines(x=c(quantile(boot.denounce$time_year,probs=.025),quantile(boot.denounce$time_year,probs=.975)),y=c(4,4))
points(mean(dat$denounce[which(dat$time==0)]),4,pch=19)
axis(2,at=c(4,5,6),labels=c("More Than A Year","One Month","Yesterday"),las=1)
# month
lines(x=c(quantile(boot.denounce$time_month,probs=.025),quantile(boot.denounce$time_month,probs=.975)),y=c(5,5))
points(mean(dat$denounce[which(dat$time==1)]),5,pch=19)
# day
lines(x=c(quantile(boot.denounce$time_now,probs=.025),quantile(boot.denounce$time_now,probs=.975)),y=c(6,6))
points(mean(dat$denounce[which(dat$time==2)]),6,pch=19)
## LOCATION
# imprecise:
lines(x=c(quantile(boot.denounce$loc_imprecise,probs=.025),quantile(boot.denounce$loc_imprecise,probs=.975)),y=c(1,1))
points(mean(dat$denounce[which(dat$loc==0)]),1,pch=19)
axis(2,at=c(1,2),labels=c("Unknown Location","Exact Location"),las=1)
# precise:
lines(x=c(quantile(boot.denounce$loc_precise,probs=.025),quantile(boot.denounce$loc_precise,probs=.975)),y=c(2,2))
points(mean(dat$denounce[which(dat$loc==1)]),2,pch=19)
box()
### Sanction
plot(c(NA,NA),xlim=c(-.1,3.1),ylim=c(.9,9.1),main="Impose Economic Sanctions",axes=F,xlab="",ylab="")
axis(1,at=c(0,1,2,3),labels=c("Strongly Oppose","Oppose","Support","Strongly Support"))
abline(v=quantile(dat$sanction,.25),lty=2,col="grey80") # quantile, all treatment groups
abline(v=quantile(dat$sanction,.75),lty=2,col="grey80") # quantile, all treatment groups
abline(v=mean(dat$sanction),lty=2) # sample mean, all treatment groups
## NUMBER
# imprecise:
lines(x=c(quantile(boot.sanction$num_imprecise,probs=.025),quantile(boot.sanction$num_imprecise,probs=.975)),y=c(8,8))
points(mean(dat$sanction[which(dat$num==0)]),8,pch=19)
axis(2,at=c(8,9),labels=c("Unknown Number","Exact Number"),las=1)
# precise:
lines(x=c(quantile(boot.sanction$num_precise,probs=.025),quantile(boot.sanction$num_precise,probs=.975)),y=c(9,9))
points(mean(dat$sanction[which(dat$num==1)]),9,pch=19)
## TIME
# year
lines(x=c(quantile(boot.sanction$time_year,probs=.025),quantile(boot.sanction$time_year,probs=.975)),y=c(4,4))
points(mean(dat$sanction[which(dat$time==0)]),4,pch=19)
axis(2,at=c(4,5,6),labels=c("More Than A Year","One Month","Yesterday"),las=1)
# month
lines(x=c(quantile(boot.sanction$time_month,probs=.025),quantile(boot.sanction$time_month,probs=.975)),y=c(5,5))
points(mean(dat$sanction[which(dat$time==1)]),5,pch=19)
# day
lines(x=c(quantile(boot.sanction$time_now,probs=.025),quantile(boot.sanction$time_now,probs=.975)),y=c(6,6))
points(mean(dat$sanction[which(dat$time==2)]),6,pch=19)
## LOCATION
# imprecise:
lines(x=c(quantile(boot.sanction$loc_imprecise,probs=.025),quantile(boot.sanction$loc_imprecise,probs=.975)),y=c(1,1))
points(mean(dat$sanction[which(dat$loc==0)]),1,pch=19)
axis(2,at=c(1,2),labels=c("Unknown Location","Exact Location"),las=1)
# precise:
lines(x=c(quantile(boot.sanction$loc_precise,probs=.025),quantile(boot.sanction$loc_precise,probs=.975)),y=c(2,2))
points(mean(dat$sanction[which(dat$loc==1)]),2,pch=19)
box()
### Intervention
plot(c(NA,NA),xlim=c(-.1,3.1),ylim=c(.9,9.1),main="Send Troops",axes=F,xlab="",ylab="")
axis(1,at=c(0,1,2,3),labels=c("Strongly Oppose","Oppose","Support","Strongly Support"))
abline(v=quantile(dat$interven,.25),lty=2,col="grey80") # quantile, all treatment groups
abline(v=quantile(dat$interven,.75),lty=2,col="grey80") # quantile, all treatment groups
abline(v=mean(dat$interven),lty=2) # sample mean, all treatment groups
## NUMBER
# imprecise:
lines(x=c(quantile(boot.interven$num_imprecise,probs=.025),quantile(boot.interven$num_imprecise,probs=.975)),y=c(8,8))
points(mean(dat$interven[which(dat$num==0)]),8,pch=19)
axis(2,at=c(8,9),labels=c("Unknown Number","Exact Number"),las=1)
# precise:
lines(x=c(quantile(boot.interven$num_precise,probs=.025),quantile(boot.interven$num_precise,probs=.975)),y=c(9,9))
points(mean(dat$interven[which(dat$num==1)]),9,pch=19)
## TIME
# year
lines(x=c(quantile(boot.interven$time_year,probs=.025),quantile(boot.interven$time_year,probs=.975)),y=c(4,4))
points(mean(dat$interven[which(dat$time==0)]),4,pch=19)
axis(2,at=c(4,5,6),labels=c("More Than A Year","One Month","Yesterday"),las=1)
# month
lines(x=c(quantile(boot.interven$time_month,probs=.025),quantile(boot.interven$time_month,probs=.975)),y=c(5,5))
points(mean(dat$interven[which(dat$time==1)]),5,pch=19)
# day
lines(x=c(quantile(boot.interven$time_now,probs=.025),quantile(boot.interven$time_now,probs=.975)),y=c(6,6))
points(mean(dat$interven[which(dat$time==2)]),6,pch=19)
## LOCATION
# imprecise:
lines(x=c(quantile(boot.interven$loc_imprecise,probs=.025),quantile(boot.interven$loc_imprecise,probs=.975)),y=c(1,1))
points(mean(dat$interven[which(dat$loc==0)]),1,pch=19)
axis(2,at=c(1,2),labels=c("Unknown Location","Exact Location"),las=1)
# precise:
lines(x=c(quantile(boot.interven$loc_precise,probs=.025),quantile(boot.interven$loc_precise,probs=.975)),y=c(2,2))
points(mean(dat$interven[which(dat$loc==1)]),2,pch=19)
legend(x="bottomright",legend=c("Group Mean (95\\% CI) \\ \\ \\ \\ ","Sample Mean","Quartiles"),
       pch=c(19,NA,NA),lty=c(1,2,2),col=c("black","black","grey80"),xpd=F)
box()
dev.off()


# end log
#sink()
